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Abstract 

We study the dynamics of network-coupled phase oscillators in the presence of coupling frustration. It was recently 
demonstrated that in heterogeneous network topologies, the presence of coupling frustration causes perfect phase 
synchronization to become unattainable even in the limit of infinite coupling strength. Here, we consider the important 
case of heterogeneous coupling functions and extend previous results by deriving analytical predictions for the total 
erosion of synchronization. Our analytical results are given in terms of basic quantities related to the network structure 
and coupling frustration. In addition to fully heterogeneous coupling, where each individual interaction is allowed 
to be distinct, we also consider partially heterogeneous coupling and homogeneous coupling in which the coupling 
functions are either unique to each oscillator or identical for all network interactions, respectively. We demonstrate the 
validity of our theory with numerical simulations of multiple network models, and highlight the interesting effects that 
various coupling choices and network models have on the total erosion of synchronization. Finally, we consider some 
special network stmctures with well-known spectral properties, which allows us to derive further analytical results. 


1. Introduction 

Self-organization and emergent collective behavior 
represent universal concepts that are vital in many non¬ 
linear processes HUE). Synchronization of large ensem¬ 
bles of coupled oscillators plays a particularly impor¬ 
tant role in our understanding of complex and network- 
coupled dynamical systems 0. Examples of the im¬ 
portance of synchronization can be found in natural 
phenomena, for instance the functionality of cardiac 
pacemakers El, mammalian circadian rhythms 0, and 
rhythmic flashing of fireflies 0, as well as engineered 
systems, for instance arrays of Josephson junctions 0, 
the power grid |8), and pedestrian bridges 0. A par¬ 
ticularly useful model for studying the synchroniza¬ 
tion of nonidentical oscillators was developed by Ku- 
ramoto ED- who showed that under suitable conditions 
the dynamics of N coupled oscillators can be reduced 
to the dynamics of N phase angles 0;, for i = 1,..., N. 
When placed on a network whose structure dictates the 
oscillators’ interaction patterns, the evolution of each 
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phase is given by 

N 

k = Wi + K ]T - 9 t ), Cl) 

i=i 

where the natural frequency uii describes the prefered 
frequency of oscillator i in the absence of coupling, 
AT > 0 is the global coupling strength, the adjacency 
matrix Aij encodes the network interactions, which is 
assumed to be undirected such that A t j = A ?i , and 
Hij ( 6 ) is the coupling function that describes the func¬ 
tional effect of oscillator j on oscillator i, which is 
assumed to be 27r-periodic and continuously differen¬ 
tiable. 


The dynamics exhibited by Eq. (|T]) have been studied 
in various contexts mm and have advanced our un¬ 
derstanding of collective behavior, particularly regard¬ 
ing the interplay between structure and dynamics and 
their effects on synchronization. Typically, the extent of 
phase synchronization of the oscillators is measured by 
the classical Kuramoto order parameter r that is defined 
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r ^ = N^ ei<>j ’ (2) 

3 = 1 

where the complex number re represents the oscil¬ 
lators’ centroid in the complex unit circle. In particu¬ 
lar, the order parameter r ranges from 0 to 1, indicat¬ 
ing complete incoherence and perfect synchronization, 
respectively, while intermediate values typically corre¬ 
spond to partial synchronization. Alternatively, sev¬ 
eral studies have defined the degree of phase synchro¬ 
nization using a combination of the collection of local 
order parameters, defined ne' l ^ i = Aije iej for 

i = 1,... ,N H3][23]|27). 

A key element of the model in Eq. (jT} is the choice of 
coupling functions fT.y(0) that defines the interactions 
between oscillators. For instance, the choice Hij(9) = 
sin (6) yields the classical Kuramoto model flOl , while 
the presence of additional modes can give rise to multi¬ 
branch entrainment, a.k.a. cluster synchronization lf30l - 
E3- Here, we focus our attention on systems with cou¬ 
pling frustration, as indicated by one or more non-zero 
values of the quantity 

h ij = H ZJ (0)/V2H' ] (0). (3) 

The physical interpretation of coupling frustration cor¬ 
responds to the case where the networks’ interaction 
terms do not all vanish when all phases are equal. The 
presence of coupling frustration is vital in the model¬ 
ing of excitable and reaction-diffusion dynamics for the 
reason that neighboring elements typically do not react 
simultaneously, but rather one after another Ii38l . Many 
such examples exist in biological and chemical systems, 
including neuron excitation ll39l . cardiac dynamics fl40l . 
and the Belousov-Zhabotinsky reaction ED- Addition¬ 
ally, coupling frustration has been linked to the emer¬ 
gence of chimera states RSlBOll . non-universal synchro¬ 
nization transitions ED, and other effects 1521 . 

In a recent publication lt53l we reported a novel phe¬ 
nomenon for networks of coupled oscillators that we 
called erosion of synchronization. In particular, we 
found that in the presence of both coupling frustration 
and structural heterogeneity the perfectly synchronized 
state (i.e., r = 1 , or equivalently, 9i = 9 2 = • • • = On) 
becomes unattainable in steady-state even in the limit of 
infinite coupling strength. To quantify the total erosion 
of synchronization in a network, we consider the quan¬ 
tity 1 — r in the limit K —> 00 , denoted 1 — r°°. We 
demonstrated this by considering the case of homoge¬ 


neous coupling, i.e., = H{9), and subsequently 

showed that the total erosion of synchronization could 
be separated into the product of two terms describing 
the contributions of coupling frustration and structural 
heterogeneity, respectively, and that both of these terms 
amplify the total erosion of synchronization. 

In this Article, we provide a more complete descrip¬ 
tion of this phenomenon. In particular, we extend our 
previous results to account for the important case of het¬ 
erogeneous coupling, i.e., when the coupling function 
governing the interaction between each pair of network 
neighbors may be distinct. We refer to this most gen¬ 
eral case, where each Hij(6) is potentially different, as 
full coupling heterogeneity . In this case we assume that 
each undirected link has an associated coupling func¬ 
tion, so that = Hji{6). We also treat the case 

where each oscillator has its own coupling function, i.e., 
Hij(6 ) = Hi(9), which we refer to as partial coupling 
heterogeneity . Unlike the homogeneous coupling case, 
in both the fully and partially heterogeneous coupling 
cases we find that the total erosion of synchronization 
cannot be separated into a product of contributions from 
the coupling frustration and structural heterogeneity. 

The remainder of this Article is organized as follows. 
In Section [2] we present our theoretical results, which 
extend previous results for homogeneous coupling to 
the cases of both full and partial coupling heterogene¬ 
ity. In Section[3]we present results from numerical sim¬ 
ulations that support our theory and explore the inter¬ 
play between coupling frustration and structural hetero¬ 
geneity. In Section [4] we study the stability of the syn¬ 
chronized state. In Section [5] we investigate erosion of 
synchronization in several network models with well- 
known spectral properties, allowing us to develop fur¬ 
ther analytical results. In particular, we consider the 
star and chain networks, as well as Watts-Strogatz net¬ 
works EH- Finally, in Section [6] we conclude with a 
discussion of our results. 

2. Theory 

In this Section we present a theoretical framework for 
quantifying the erosion of synchronization for the dy¬ 
namics defined in Eq. 0- We begin by considering the 
case of fully heterogeneous coupling, i.e., where each 
undirected link connecting oscillators i and j have a po¬ 
tentially distinct coupling function H,,j(9) = Hji(9). 
We also consider the case of partially heterogeneous 
coupling, i.e., when each oscillator has its own coupling 
functions, Hij(6) = Hi{9). Finally, we compare these 
results to the originally derived results for homogeneous 
coupling, i.e., 11^(9) = H(9 ), presented in Ref. {53l . 
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2.1. Fully heterogeneous coupling 

We begin by following Ref. l26l and consider the dy¬ 
namics of Eq. 0 in the strong coupling regime, i.e., 
r ss 1. In typical networks, such a state can be attained 
in a variety of ways, most readily by considering either 
a sufficiently large coupling strength, or a set of nat¬ 
ural frequencies with a sufficiently small spread. It is 
worth pointing out that these two situations are equiv¬ 
alent up to a rescaling of time, and thus the results 
presented here are valid in both cases. In the strong 
coupling regime the oscillators become tightly packed 
around the mean phase i/j, implying that 10,; — 6j <C 1 
for all (i,j) pairs. Thus, the contribution of each pair¬ 
wise interaction can be linearized to Hij{9j — 9i) k, 
Hij( 0) -f H l i j(0)(9j — 9i), and Eq. 0 can be approxi¬ 
mated by 


N 

9/ ~ tUj -f-A d'i : A ^ ^ Lij9j • 
i=i 


or rather in vector form, 

0 K,oj + I<d- KLG. 


(4) 


(5) 


Here, d and L represent the weighted degree vector 
and weighted Laplacian matrix. In contrast to the un¬ 
weighted degree vector d and unweight Laplacian ma¬ 
trix L , whose entries are defined 

N 

d-i — ^ -4 fj. Lij — Stjdi ij , ( 6 ) 

i=i 

the entries of their weighted counterparts are given by 


N 


dj — ^ ^ F,j (0), 


(7) 


i=i 


N 


Lij = Sij E (°) ^ A ^ H ij (°). (8) 


i=l 


where Sij is the Kronecker delta. 

We now aim to solve for the steady-state of Eq. 0 
where, in a synchronized state, each oscillator trav¬ 
els at the same effective frequency, (I, = 51 for all 
i = 1 ,... ,N. This effective frequency, which is also 
given by FI = (9) = N _1 JA 9i is simply the mean of 
the heterogeneous part of Eq. 0, i.e., FI = (w + Kd) = 
N- 1 + Kdi, which follows from the fact that 

L is symmetric and maps vectors to the space of vec¬ 
tors with zero mean (see explanation below). Insert¬ 
ing 0 = Vl = fl[l,..., 1] T in Eq. (j5j and rearrang¬ 


ing, we effectively enter the rotating reference frame 
9 h -> 6 + Fit and obtain 

LG* = (u> + Kd - FI)/K, (9) 


where G* denotes the steady-state solution in the rotat¬ 
ing frame. 

To solve for G*, we now aim to define the pseudoin¬ 
verse of the weighted Laplacian matrix ea, which can 
be defined using the spectral properties of L. Assum¬ 
ing that the network is connected and undirected, L has 
a single zero eigenvalue with the remainder being real 
and positive so that they can be ordered 0 = Ai < 
A 2 < • • • < Ay. Furthermore, the eigenvector that cor¬ 
responds to the trivial eigenvalue Ai = 0 is simply the 
constant vector v l oc [1,..., 1] T and represents motion 
along the synchronization manifold. With the remainder 
of the eigenvectors (which we assume are normalized to 
||i> J '|| = 1) the pseudoinverse can be defined 


N 


zt = E 

3 =2 


V^V^ T 

— 


( 10 ) 


Importantly, we note that L and lJ share the nullspace 
of the all vectors spanned by v 1 , and thus map all vec¬ 
tors to the N — 1-dimensional space of vectors in R N 
with zero mean. 

We can now return to Eq. (|9]) and apply the pseudoin¬ 
verse to obtain 


e* = D y , (ii) 


where 


y = u;/K + d. 


( 12 ) 


We note that $7 does not appear in Eq. (12 1 since II 
0. To evaluate the degree of synchronization, we now 
consider the order parameter given by Eq. 0- We note 
that by a suitable shift in initial conditions, the mean 
phase %l> can be set to zero without loss of generality. 
Noting that both sides of Eq. 0 must be purely real, 
the right-hand side of Eq. 0 can be expanded to obtain 


3* II2 


r ~ 1 — 


2 TV 


(13) 


Inserting Eq. 0 into Eq. 0, noting that \\G* || 2 = 
(G*, G*), and using the definition of L f in Eq. (101, we 
obtain 


r ~ 1 — 


J(y,L) 

2 


(14) 
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where J is the synchrony alignment function previously 
studied in Refs. 12611531 and is defined by 




if W 

N h 


(15) 


Finally, to evaluate the total erosion of synchronization 
we consider the deviation from perfect synchronization, 
1 — r, in the limit K —> oo, which we denote 1 — r°°. 
In this limit we have that y = d, and we finally obtain 


J(d, L) 

'~2 


(16) 


Equation ( fl6| quantifies the total erosion of synchro¬ 
nization in a given network. We point out that the con¬ 
tributions from the coupling frustration [i.e., Hjj (0) and 
H' J (0)] and the network structure (i.e., d and L) do not 
separate, as they jointly define the entries of both the 
weighted degree vector d and the weighted Laplacian 
L. 


2.2. Partially heterogeneous coupling and homoge¬ 
neous coupling 


Next we consider two special cases of coupling: par¬ 
tially heterogeneous coupling and homogeneous cou¬ 
pling. In the case of partial heterogeneity we assume 
that each oscillator has its own coupling function, i.e., 
Hij{6) = Hi{9). Carrying through an analysis similar 
to that outlined in Eqs. (|4|-([T6|, we find that 


1 - 


J(d,L) 

'~2 


(17) 


where 


Table 1: Summary of coupling types, total erosion of synchronization, 
and separability. Definition of the 
synchrony alignment function J is given in Eq. (fl5j). 


Coupling Type 

1 — r°° 

Horn., Hijf) = Hf) 

Par. Het., iT.y-(-) = H t {-) 
Fully Het., H vj ( : ) 

2 jjJ(o) L ) (separable) 

J(d, L)/2 (not separable) 
J{d, L)/2 (not separable) 


Here, the right-hand side of Eq. (19 1 separates conve¬ 
niently into the product of the square of the coupling 
frustration h = \H (ft) / \/2H'{0)\ and structural hetero¬ 
geneity in the network, as measured by J(d,. L). This 
separation is a convenient property that is reminiscent 
of the separation of dynamics and structure in the anal¬ 
ysis of synchronization of identical and nearly identi¬ 
cal oscillators using the Master Stability Function ap¬ 
proach I56l,l57) . 

For convenience, we summarize in Table |T| the ex¬ 
pression for the total erosion of synchronization for 
each type of coupling we have considered, and indicate 
whether or not the contribution from coupling frustra¬ 
tion and network structure are separable. We also note 
that in each of the three types of coupling we have con¬ 
sidered here, we find that if no coupling frustration is 
present, i.e., all // )y (0) = 0, then the total erosion of 
synchronization is zero, i.e., 1 — r°° = 0. This fol¬ 
lows directly in the homogeneous case [Eq. ( [T9| ] from 
the coefficient H (0). In the fully and partially heteroge¬ 
neous cases [Eqs. (j9j) and ( [To| ] //,, (0) = 0 implies that 
d = 0 and d = 0, respectively, and in turn J(0, L) = 0 
for any choice of L. 


3. Numerics 


di = Hi^di/H'iO). (18) 

Similar to the case of full heterogeneity, the contribu¬ 
tions to the total erosion of synchronization from cou¬ 
pling frustration and network structure do not sepa¬ 
rate, this time due only to the construction of the new 
weighted degree vector d. 

Finally, in the case of homogeneous coupling, each 
coupling function is identical, i.e., Hij(9) — H(6). 
This was the case studied in Ref. lf53l , and the theory 
outlined in Eqs. 0-<[T6]) simplifies even more so, yield¬ 
ing 


1 — r 


OO 


7T 2 (0) 

2H' 2 (0) 


J(d,L). 


(19) 


In this Section we illustrate the essential properties 
of the theory outlined above using several numerical 
examples. As a benchmark choice of coupling, from 
this point forward we consider Sakaguchi-Kuramoto- 
type coupling given by a sinusoid with an associated 
phase-lag, i.e., Hij{6 ) = sin(0 — a*j) l58l . In par¬ 
ticular, the parameter represents a phase-lag be¬ 
tween oscillators i and j , and can be taken from the 
range a,,- £ (—7 t/2,7t/ 2) (although here we will re¬ 
strict our focus on only non-negative values). In the 
case of full heterogeneity, each undirected link (i. j) is 
then assigned its own phase-lag a tJ = aji, whereas in 
the case of partial heterogeneity each oscillator i is as¬ 
signed its own phase-lag a r . In the homogeneous case 
where a single phase-lag is chosen for all interactions. 
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Figure 1: (Color online) Erosion of synchronization: Effect of coupling heterogeneity and frustration. Deviation from perfect synchronization 1 — r 
vs coupling strength K in a SF network of size N = 1000, exponent 7 = 3, and mean degree ( d) = 4. Results are presented for (a) homogeneous 
coupling, (b) partially heterogeneous coupling, and (c) fully heterogeneous coupling for phase-lag parameter a = 0 (black circles), 0.05 (blue 
triangles), 0.1 (green squares), 0.15 (orange inverted triangles), and 0.2 (red diamonds). Results for partially and fully heterogeneous couplings 
are averaged over 10 independent phase-lag realizations. The horizontal dashed lines in the right, center and left columns indicate the theoretical 
predictions for 1 — r°° given by Eqs. {76} , fT?} and {79}. respectively. 


the coupling frustration is described by a single param¬ 
eter h = \H(0)/^/2H'(0)\ = | tan(— a)\/s/2. 

For all simulations presented in the remainder of this 
paper, we solve Eq. 0 using a two-step Runge-Kutta 
method, a.k.a., Heun’s method l59l . We use a time 
step At = 4 x 1CT 4 , maintaining numerical stability 
for large coupling strengths by rescaling time by Ff _1 . 
Steady-state is obtained by discarding a significant tran¬ 
sient of at least 10 5 time steps. In each simulation we 
draw natural frequencies independently from the unit 
normal distribution. We will now present numerical re¬ 
sults illustrating the effects of different phase-lags and 
network structures. 

3.1. Effect of phase-lags 

We begin by considering the effect of different phase- 
lags on the erosion of synchronization for a given net¬ 
work. We first consider scale-free networks with power- 
law degree distribution P(d) ex g? - 7 for d > cLq built 
using the configuration model |[60l . In order to tune the 
mean degree (d) of each network, we set the minimum 
degree equal to do = (7 — 2)(d)/(7 — 1). We assign 
phase-lags according to a mean phase-lag parameter a 
as follows. In the case of homogeneous coupling, a 
describes the global phase-lag for the entire network. 
Otherwise, each distinct phase-lag {at or a l:l ) is drawn 
independently and uniformly from the interval [0, 2a] to 
maintain a mean phase-lag of a. 

We begin our investigations by considering a SF net¬ 
works size N = 1000 with exponent 7 = 3 and 
mean degree (d) = 4, and simulating the dynamics 
of Eq. 0 for homogeneous, partially heterogeneous, 
and fully heterogeneous coupling with several a values. 
We plot the resulting profiles 1 — r vs K for each type 


of coupling in Figs. [I] a), (b), and (c), respectively, us¬ 
ing a log-log scale. Results using phase-lag parameters 
a = 0, 0.05, 0.1, 0.15, and 0.2 are plotted using black 
circles, blue triangle, green squares, orange inverted tri¬ 
angles, and red diamonds, respectively. For the partially 
and fully heterogeneous cases, results represent an av¬ 
erage over 10 simulations with different phase-lag real¬ 
izations. First, when no coupling frustration is present 
in the system (i.e., a = 0), the quantity 1 — r decays as a 
power-law as K is increased in each panel. This decay 
is observed to continue well past the windows shown. 
For nonzero values of a, each curve 1 — r initially de¬ 
cays, then saturates at a finite positive value, indicating 
that the perfectly synchronized state r = 1 cannot be 
attained. For each of the three types of coupling hetero¬ 
geneity, larger phase-lags yield a greater total erosion 
of synchronization. For each nonzero value of a, we 
also plot the theoretical saturating value of 1 — r°° in 
horizontal dashed lines [as given by Eqs. ( fl6] >, 07} anc * 
©], which accurately predict the results from simula¬ 
tions. We also observe that the total erosion of synchro¬ 
nization tends to be larger on aggregate for partially het¬ 
erogeneous coupling than for either fully heterogeneous 
or homogeneous coupling. 


3.2. Effect of network structure 

Next we consider the effects of network structure on 
the erosion of synchronization. Using the same net¬ 
work model as in Sec. 3.1 we first explore different 
topologies by varying the power-law exponent 7 . In 
Figs. 0a), (b), and (c) we plot 1 — r vs K for networks 
of size N = 1000 with mean degree (d) = 4 for ho¬ 
mogeneous, partially heterogeneous, and fully hetero¬ 
geneous coupling, respectively, for fixed a = 0.2. To 
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Figure 2: (Color online) Erosion of synchronization: Effect of network structure. Deviation from perfect synchronization 1 — r vs coupling strength 
K in various networks of size N = 1000 and mean degree (d) = 4 using a fixed phase-lag parameter a = 0.2. Panels (a)-(c) display results 
for SF networks with homogeneous, partially heterogeneous, and fully heterogeneous coupling, respectively, and exponent 7 = 4.5 (blue circles), 
4 (green triangles), 3.5 (orange squares), and 3 (red diamonds). Panels (d)-(f) display results for BAER networks with homogeneous, partially 
heterogeneous, and fully heterogeneous coupling, respectively, and heterogeneity parameter (3 = 0 (blue circles), 0.33 (green triangles), 0.67 
(orange squares), and 1 (red diamonds). Each data point represents an average over 10 networks and 10 independent phase-lag realizations for 
partially and fully heterogeneous couplings. The horizontal dashed lines in the right, center and left columns indicate the theoretical predictions 
for 1 — r°° given by Eqs. m , ® and 0 - respectively. 


vary the degree heterogeneity in each network, we use 
networks with corresponding exponent 7 = 4.5, 4, 3.5, 
and 3, the results for which are plotted in blue circles, 
green triangles, orange squares, and red diamonds, re¬ 
spectively. Each curve represents the average over 10 
network realizations and phase-lag draws. For this net¬ 
work model we observe that networks with larger de¬ 
gree heterogeneity (smaller 7 ) yield larger total erosion 
of synchronization for all three types of coupling. 

We contrast these results with results for another 
network model that allows for tunable degree hetero¬ 
geneity. Specifically, we consider a network growth 
model ISTI that interpolates between Erdos-Renyi (ER) 
ll62l type random networks and Barabasi-Albert (BA) 
l63l preferential attachment type networks, which we 
summarize as follows. We prescribe a heterogeneity pa¬ 
rameter jd £ [ 0 , 1 ] and a minimum degree do, and start 
with a fully-connected set of do 4- 1 nodes. Next, nodes 
are added one-by-one to the network, each making do 
links to the existing nodes. Each link can be made in one 
of two ways: with probability (3 the link is made prefer¬ 
entially such that the new node connects to an existing 
node i with probability pi oc di — to, and otherwise the 
link is made to a node uniformly at random. The param¬ 
eter to < do affects the final degree distribution. Nodes 


are added to the network until the desired size of the net¬ 
work N is attained. If j3 = 0, all links are constructed 
uniformly at random, and the resulting network has an 
ER-like topology. If jd = 1, the process is equivalent to 
a BA preferential attachment method, resulting in an SF 
structure with degree distribution 7 « 3 — mjdo- Inter¬ 
mediate values of /3 result in topologies that interpolate 
these two extreme cases. We refer to this model as the 
BAER model from this point forward. Importantly, the 
heterogeneity of the resulting networks can be tuned us¬ 
ing the parameter fd, with large (small) values yielding 
more (less) heterogeneous degree distributions. 

We now explore the effect that tuning the hetero¬ 
geneity parameter fd in the BAER model has on ero¬ 
sion of synchronization. In Fig. |2fd), (e), and (f), we 
plot 1 — r vs K for BAER networks with homoge¬ 
neous, partially heterogeneous, and fully heterogeneous 
coupling, respectively. Results are shown with fixed 
a = 0.2 for networks of size N = 1000 with mean 
degree (d) = 4 (do = 2) and m = 1.6. We vary the 
heterogeneity parameter in each network, plotting the 
results from /3 = 0, 0.33, 0.67, and 1 in blue circles, 
green triangles, orange squares, and red diamonds, re¬ 
spectively. Each curve represents the average over 10 
network realizations and phase-lag draws. In contrast 
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Figure 3: (Color online) SF vs BAER network topologies. Theoretically predicted total erosion of synchronization 1 — r°° vs the standard deviation 
of the degree distribution a = yj{d 2 ) — ( d ) 2 for SF (blue circles) and BAER (red triangles) networks of size N = 500 with mean degree 
(d) = 4 using a phase-lag parameter a = 0.2. Panels (a)-(c) display results for homogeneous, partially heterogeneous, and fully heterogeneous 
couplings, respectively. Each data point represents an average over 1000 network realizations and 10 independent phase-lag realizations for 
partially and fully heterogeneous couplings. 


to the results obtained with SF networks (see the top 
row of Fig. [2] wherein increased heterogeneity yielded 
an increase in the total erosion of synchronization for 
all three coupling heterogeneities), we observe the op¬ 
posite effect for BAER networks. In particular, as de¬ 
gree heterogeneity is increased (i.e., 13 is increased) the 
total erosion of synchronization decreases for all three 
types of coupling heterogeneity. Finally, we note that 
similar to our observation made for Fig. □ in Fig. 0 one 
can also observe that the total erosion of synchroniza¬ 
tion tends to be larger in the case of partially hetero¬ 
geneous coupling than for either fully heterogeneous or 
homogeneous coupling. 

We next explore more closely the relationship be¬ 
tween total erosion of synchronization and degree het¬ 
erogeneity in the SF and BAER network models. 
Specifically, we compute for an ensemble of networks 
the predicted value 1 — r°° and the standard deviation of 
the degree distribution a = i J ( d 2 ) — (d) 2 for various 7 
and p and fixed (d) = 4 and a = 0.2. In Fig.[3]we plot 
1 — r°° vs a for SF and BAER networks (blue circles 
and red triangles, respectively). Results are shown for 
homogeneous, partially heterogeneous, and fully het¬ 
erogeneous coupling in panels (a), (b), and (c), respec¬ 
tively. Each data point represents the average over 1000 
networks of size N = 500, each with 10 independent 
phase-lag realizations for the cases with heterogeneous 
coupling. In all three panels, as er increases, the total 
erosion of synchronization increases for the SF model 
but decreases for the BAER model. 

4. Stability 

In this Section we shift our focus to the stability of 
the stationary state given in Eq. ( fTT) . We consider here 
the limit of strong coupling strength, i.e., K -A 00 . The 


linear stability of this state is dictated by the spectrum 
of the scaled Jacobian matrix DF = DF/K , whose 
entries are given by 

nr - J A l2 H[ 2 (Q)~QX) if i + j 

lJ l - Ej/i AjH'v (e* 0 :) if * = j. 

( 20 ) 

Like the Laplacian matrix, each row of DF sums to 
zero, yielding a trivial zero eigenvalue that corresponds 
to the translational invariance of the dynamics in Eq. Q- 
A stationary solution 6* is then linearly stable if the real 
parts of all nontrivial eigenvalues are negative, or con¬ 
versely, if no eigenvalue has positive real part. More 
precisely, given the eigenvalues A i of DF we define 
A™“ = maxj{Re(Aj) : A, ^ 0}, such that A™ x < 0 
(> 0 ) implies that the state 6 * is linearly stable (unsta¬ 
ble). 

To guide our analysis, we first consider the case of 
no coupling frustration, 7/^(0) = 0 , for which 9* = 0 . 
This follows from the fact that each entry of the vector d 
is zero and implies that the scaled Jacobian is negatively 
proportional to the unweighted Laplacian, DF oc —L. 
Assuming that the network is strongly connected, the 
nontrivial eigenvalues of L all have positive real part, 
implying that the nontrivial eigenvalues of DF have all 
negative real part, and thus the solution is stable. 

Next, we consider solutions under the presence of 
coupling frustration. Given that there exists a spectral 
gap for the unweighted Laplacian, a sufficiently small 
amount of coupling frustration can be found such that 
the solution 8* is guaranteed to remain stable. The 
extent to which coupling frustrations can be increased 
without destabilizing the solution, however, is unclear. 
To shed some light on this, we consider the structure 
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of DF. For the typical case of // (; (0) > 0 and small 
enough 1 9* — 0*\, the scaled Jacobian consists of nega¬ 
tive diagonal entries with non-negative off-diagonal en¬ 
tries. In particular, as long as this structure is main¬ 
tained, DF remains negative semi-definite and its spec¬ 
trum remains bounded in the left-half complex plane. 
(This can be shown, for instance, by using the Gersh- 
gorin circle theorem ED .) A necessary condition for 
instability is then given by at least one off-diagonal en¬ 
try becoming negative, which allows for the possibil¬ 
ity of one or more eigenvalues crossing the imaginary 
axis. In the case of homogeneous coupling, we showed 
in Ref. Il53l that a necessary condition for the loss of 
linear stability is given by 

(21) 


For partially and fully heterogeneous coupling, this con¬ 
dition becomes, respectively, 

min H[ (^[L^d]j — [L^d]^ < 0, (22) 


and 


juiin H-j ([Dd]j — [Dd]^ < 0. (23) 

When the appropriate equation of Eqs. ([2T|)-([23| is false 
the synchronized state is guaranteed to be linear stable, 
but if it holds true, this admits the possibility that it is 
unstable. 

To illustrate the utility of Eqs. (|2T|-([23 |i, we study the 
stability BAER networks of size N = 500 with /3 = 0 
with homogeneous, partially heterogeneous, and fully 
heterogeneous coupling. (We note that results obtained 
using other parameters and network models were found 
to be similar.) First, using a sample of 1000 networks 
each with 10 independent realizations of phase-lags for 
each value of a, we plot in Fig. |4|a) the fraction of in¬ 
stances where the synchronized state was found to be 
unstable. Results for homogeneous, partially heteroge¬ 
neous, and fully heterogeneous coupling are plotted us¬ 
ing blue circles, red triangles, and green squares, respec¬ 
tively. We also calculate the critical value of a for which 


the conditions given in Eqs. (21 1 —(23 i first become true 


on average, which we denote with using thick vertical 
markers. Our results indicate that these conditions do 
in fact provide a reasonably good lower bound for the 
transition to instability. 

In Fig. |4jb) we plot the average real part of the 
largest nontrivial eigenvalue A™ ax of DF calculated 



Figure 4: (Color online) Stability of the synchronized solution. For 
BAER networks of size N = 500 and 0 = 0, (a) the fraction of so¬ 
lutions found to be unstable for homogeneous (blue circles), partially 
heterogeneous (red triangles), and fully heterogeneous (green squares) 
coupling. Thick vertical markers indicate the critical value a where 
the necessary instability conditions given by Eqs. iJTTJ— (??). All re¬ 
sults indicate mean values computed across 100 simulations per each 
a value, (b) The mean real part of the largest nontrivial eigenvalue 
(max A re ) (symbols) plus and minus (dashed curves) the standard de¬ 
viation. Results represent the average over 100 network realizations 
and 10 independent phase-lag realizations for partially and fully het¬ 
erogeneous coupling. 


from these simulations, denoting the standard devia¬ 
tions using dashed curves. Recall that for a given sys¬ 
tem, A™“ < 0 and A“ ax > 0 imply stability and instabil¬ 
ity, respectively. Therefore, the transitions from stabil¬ 
ity to instability shown in panel (a) are driven by the de¬ 
pendence of A™ ax on a. Naturally, for small a the values 
of A™ x do not depend strongly on the type of coupling 
heterogeneity. However, as a increases, the dependence 
of A™ ax on a is varies for each type of coupling hetero¬ 
geneity, which leads to the different transitions shown 
in panel (a). For example, A™ ax grows faster for the par¬ 
tially and fully heterogeneous couplings than for homo¬ 
geneous coupling, which leads to a transition from sta¬ 
bility to instability that occurs for smaller a. Moreover, 
the standard deviation in A™ ax is much smaller for ho¬ 
mogeneous coupling than for the other couplings, and 
thus the transition from stability to instability is much 
more abrupt for homogeneous coupling than for the het¬ 
erogeneous couplings. 


5. Special network structures 

In all of the results presented above, both for pre¬ 
dicting the total erosion of synchronization and classi¬ 
fying the stability of the solution, a key ingredient is 
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(a) 




Figure 5: (Color online) Special network structures. Example illus¬ 
trations of (a) star, (b) chain, and (c) Watts-Strogatz networks. 


the (weighted or unweighted) network Laplacian, and 
especially its spectral properties. While the spectral 
properties for realizations of certain networks models 
can vary significantly (64), for a handful of special net¬ 
work structures the spectral properties are well-known. 
Such network structures, such as stars and chains, are 
used in various studies as either small motifs or mod¬ 
ules l l65l [661 . collections of which make up a full net¬ 
work, or on their own, serving as examples for various 
dynamical phenomena ifTTl 167 . 681. In order to obtain 
analytical predictions for the total erosion of synchro¬ 
nization in such networks, we focus here on the case 
of homogeneous coupling where the unweighted Lapla¬ 
cian and degree vector can be used. Furthermore, while 
we have observed previously that results for homo¬ 
geneous coupling differ from heterogeneous coupling 
quantitatively, we note that homogeneous coupling does 
in fact serve as a predictive benchmark, displaying qual¬ 
itatively similar results to heterogeneous coupling. 

We will begin by presenting analytical results for the 
star and chain networks, where the total erosion of syn¬ 
chronization can be predicted in terms of the network 
size. We note that all networks with regular structure, 
i.e., all nodes having the same degree, such as the ring 
or periodic lattices, trivially yield zero total erosion of 
synchronization, and therefore we forgo any more con¬ 
sideration of these networks here. Finally, we provide 
analytical approximations for the total erosion of syn¬ 
chronization in Watts-Strogatz (WS) networks (54) . an 
important network model where progress is currently 
being made for the analytical approximation of spectral 
properties (69). 


The network Laplacian matrix is then given by 



'N — 1 

-1 

-1 ... 

- 1 ' 


-1 

1 

0 ... 

0 

L = 

-1 

0 

1 ... 

0 


-1 

0 

0 ... 

1 


and the degree vector is given by d = [A r 
1,1,..., 1] T . It is easy to show that the eigenvalues 
of L are A i = 0, A 2 = = Xn-i = 1, and 

Ajy = N. The normalized eigenvectors corresponding 
to Ai and Ajv are v 1 = [1,..., 1 } T /\fN and v N = 
[TV — 1, —1,..., — l] T /v / N 2 — N and the degree vec¬ 
tor can be written as the linear combination 

d= 2N ~ 2 v 1 + (N-2)^1-N- 1 v n . (25) 

V N 

Finally, when evaluating J id,. L), only the v N part con¬ 
tributes, and after simplification we obtain 


J{d, L) 


(N — 2) 2 (1 — TV -1 ) 


(26) 


Interestingly, the behavior of J(d, L) for the star net¬ 
work is not monotonic with the network size N. At 
N = 2, J(d, L) is zero, corresponding to fact that 
the degree vector is constant, immediately above which 
J(d,L) increases. A maximum is reached at N = 
4 + 2V2, after which .1 id,. L) decreases for all larger 
N, scaling like Jid,. L) ~ N~ x . Restricting N to in¬ 
tegers greater than or equal to 2 , this implies that the 
maximum total erosion of synchronization for homoge¬ 
neous coupling in a star network is obtained at a size of 
N = 7 (i.e., six leaf nodes), above and below which the 
total erosion is less. 

In Fig. [ 6 (a) we plot our theory vs simulations us¬ 
ing star networks of various sizes and setting a = 0 . 1 . 
The theoretical curve is obtained using Eq. (26) with 
Eq. (T9) and plotted as the dashed red curve, and re¬ 
sults from simulations are plotted using blue circles. 
We note excellent agreement between theory and sim¬ 
ulation, both of which capture the non-monotonicity of 
1 — r°° depending on network size. 


5.1. Star network 

We begin by considering the star network-a network 
consisting of a single “hub” node connected to many 
“leaves” [see Fig. [5(a)]. For a star network of size N, 
we index the hub i = 1 and the leaves i = 2,..., N. 


5.2. Chain network 

Next, we consider the undirected chain-a network 
consisting of sequentially linked nodes with end nodes 
indexed i = 1 and N [see Fig.[5)b)]. The Laplacian of 
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Figure 6: (Color online) Erosion of synchronization: stars and chains. 
The total erosion of synchronization 1 — r°° as a function of network 
size N for (a) star and (b) chain networks using a = 0.1 and 0.04, re¬ 
spectively. Theoretical results given by Eq. G2> and either (a) Eq. {26| 
or (b) Eq. {29) are plotted as dashed red curves, and results from sim¬ 
ulations are plotted with blue circles. 


The dashed red curve indicates our theory, which is ob¬ 
tained using Eq. ( |29) with Eq. ( | 1 9) . and results from 
simulations are plotted using blue circles. These are in 
excellent agreement. We find that unlike star networks, 
the total erosion of synchronization for chains increases 
monotonically with increasing size N. The increasing 
behavior of J(d, L) for the chain is not only surprising 
in light of the results from the star network, but also in 
comparison to the ring. In particular, rings and chains 
of the same size differ by only the addition/subtraction 
of a single undirected link, and they in fact share similar 
spectral properties. However, whereas the total erosion 
of synchronization is trivially zero in the ring since the 
degree vector is constant, the chain yields a large total 
erosion of synchronization that increases with size. It is 
also worth noting that on the chain, the dynamics take 
a long time to relax to steady-state, which can be ex¬ 
plained by the fact that the time-scale for relaxation r is 
dictated by the inverse of the smallest non-trivial eigen¬ 
value S 2 = 4 sin 4 (n/2N), which approaches zeros as 
N approaches infinity. 


the chain is given by. 


5.3. Watts-Strogatz networks 


' 1 

-1 ... 

0 

0 ' 

-1 

2 ... 

• o 

• O 

• • o 

0 ... 

2 

-1 

1 - 

o 

0 ... 

-1 

1 _ 


and the degree vector is given by d = [1, 2 ,..., 2, l] 7 . 
It is straight forward to show that the eigenvalues of L 
are given by A j = 4sin 4 [7r(j — l)/2 N], and the cor¬ 
responding eigenvectors are v 1 = [1,..., 1 ] T /y/N and 
vi with entries given by u 4 = sj2/N cos[7r(j — l)(2z — 
1)/2IV] for j > 2. The degree vector can thus be written 
as 


d = 2 VNv 1 — e 1 — e N , (28) 


where e? is the canonical basis vector with entries e\ = 
Sij. Again, the eigenvector v 1 does not contribute to 
J(d. L), and after some algebra, we obtain 


J{d, L) 


AT 9 

i ^ ( cos Vj + cos YpA 2-^v — 1)]) 

(V2 ^ 8 sin 4 tpj ’ 

(29) 


where ipj = n(j — l)/2 N. 

In Fig. [6ja) we study the total erosion of synchro¬ 
nization for chains of various sizes using a = 0.04. 


Finally, we consider the Watts-Strogatz (WS) net¬ 
work model (54). This popular model provides an inter¬ 
polation between regular periodic lattice-type networks 
and random networks which display the small-world 
property. Here we consider the simple case of ring-like 
WS networks [see Fig. [5jc)] . Given a set of N nodes ar¬ 
ranged in a ring and a chosen uniform degree d (assum¬ 
ing to be even), each node is connected to the d/2 clos¬ 
est neighbors on each side. Next, given a rewiring prob¬ 
ability q , each link is either rewired or not rewired, with 
probabilities q and 1 — q, respectively. Each rewired link 
is then replaced with a link connecting one of the orig¬ 
inal nodes (chosen randomly) and another node that is 
chosen uniformly at random from the remaining nodes. 
Thus, for extreme values q = 0 and q = 1, the resulting 
network is a perfect ring or an Erdos-Renyi type net¬ 
work (62). 

While precise expressions for the total erosion of 
synchronization in general WS networks are difficult 
to obtain, they can be found for the limiting values of 
<7 = 0 and 1. For q = 0 the network is a regular one¬ 
dimensional lattice, and thus J (d, L) = 0 trivially since 
the degree vector d is constant. For q = 1 on the other 
hand, J{d, L) can be approximated based a simplifying 
assumption. In particular, given the completely random 
network structure, we assume that the eigenvalues A j 
and the projections of the eigenvectors onto the degree 
vector (vi, d) are uncorrelated, allowing us to separate 
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Eq. (fl9|) when q = 1. 



rewire prob., q 


Figure 7: (Color online) Erosion of synchronization: Watts-Strogatz 
networks. The total erosion of synchronization 1 — r°° as a function 
of the rewire probability q for WS networks of size N = 500 with 
mean degree ( d) = 40 using a homogeneous phase-lag a = 0.4. 
Theory given by Eq. |20| and simulations of the dynamics are plotted 
by the red dashed red curve and blue circles, respectively. We also plot 
Eq. J30J by the horizontal dashed line, which is in excellent agreement 
when q = 1. 


the synchrony alignment function: 


T(H n if («M) ! 


3 = 2 





y]«<2) 2 


(30) 


where we have used that the second term in the prod¬ 
uct is equal to the variance of the degree distribution. In 


principle, Eq. (30 1 can be used to approximate .ltd. L) 
directly from a given network’s properties (i.e., the 
degree vector and eigenvalues), or rom using known 
asymptotic results for the class of networks. For ex¬ 
ample, it is well known for ER networks that the degree 
distribution is binomial and the eigenvalue distribution 
is in certain cases well-approximated based on a semi¬ 
circle law. 

We test this approximation and investigate the behav¬ 
ior of WS networks for intermediate rewiring probabil¬ 
ities by comparing the predicted total erosion of syn¬ 
chronization given by Eq. ( |T9| ) and the results from sim¬ 
ulations for several WS networks. In Fig. [7] we plot the 
total erosion of synchronization 1 — r°° vs the rewire 
probability q for WS networks of size N = 500 with 
mean degree (d) = 40 using a = 0.4. Theoretical pre¬ 
dictions are plotted by the dashed red curve, and results 
from simulations are plotted in blue circles. Each data 
point represent an average over 20 network realizations. 
In addition, we plot by the horizontal dashed line our 
approximation for q = 1 given by Eq. which is 
in excellent agreement with both our simulations and 


6. Discussion 

In this Article we have studied erosion of synchro¬ 
nization in networks of coupled oscillators, whereby 
perfect synchronization is unattainable even in the limit 
of infinite coupling strength and is a phenomenon that 
arises in the presence of both coupling frustration and 
structural heterogeneity. We have generalized previous 
results to the important case of heterogeneous coupling, 
allowing for the interactions between different pairs of 
network neighbors to be described by different func¬ 
tions. As compared to homogeneous coupling, where a 
single coupling function describes all of the interaction 
in the network, the theoretical predictions for heteroge¬ 
neous coupling become more complicated. While the 
theoretical prediction for the total erosion of synchro¬ 
nization separates into the product of terms describing 
the coupling frustration and network structure in the ho¬ 
mogeneous case [see Eq. ©], it does not when cou¬ 
pling is heterogeneous [see Eqs. ( p~6] ) and (| 1 7[i|. How¬ 
ever, our predictions show that the presence of hetero¬ 
geneity in coupling frustrations amplifies the total ero¬ 
sion of synchronization; even when the mean frustration 
is the same as the homogeneous case, heterogeneity in 
the frustration increases the deviation from the perfectly 
synchronized state, as measured by 1 — r. Additionally, 
the heterogeneity in the coupling functions reduces the 
range of coupling frustrations for which the synchro¬ 
nized solution remains stable [see Fig. [4]. 

We have also studied erosion of synchronization in 
some special network structures. In the case of homo¬ 
geneous coupling, we have derived analytical results for 
the star and chain networks. Remarkably, in the star net¬ 
work the total erosion of synchronization is maximized 
for a star of size N = 7 (i.e., 6 leaves) and decays when 
N is either increased and decreased. In particular, for 
large N we find that 1 — r°° decays as IV -1 . In the 
case of the chain network, the total erosion of synchro¬ 
nization increases as the chain is lengthened—that is, 
1 — r°° increases with N for N < 100 (and it will have 
to decrease at some point because r°° > 0). Finally, 
we investigated the case of Watts-Strogatz networks and 
provided further analytical results for the limiting cases 
in which the rewiring probability is q = 0 or 1. 
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